#ifndef PARTICLE_SYSTEM
#define PARTICLE_SYSTEM

#include "btBulletDynamicsCommon.h"
#include "Magnet.h"
#include "Vector3.h"
#include <vector>
#include <map>
#include <math.h>
#include <algorithm>

using namespace std;

class MultiMap {

private:
	map<int, vector<int>> entries;

public:

	~MultiMap() { Clear(); };

	void Clear() { 
		//BT_PROFILE("MultiMap::Clear");
		entries.clear(); 
	};
	bool isInserted(int key) {
		//BT_PROFILE("MultiMap::isInserted");
		if (key < 0) return false;
		return entries.find(key) != entries.end();
	};
	bool Add(int key, int value) {
		//BT_PROFILE("MultiMap::Add");
		if (key < 0) return false;
		
		auto it = entries.find(key);
		if (it != entries.end()) {
			it->second.push_back(value);
			return true;
		}
		else {
			vector<int> temp;
			temp.push_back(value);
			entries.insert(pair<int , vector<int>>(key,temp));
			return true;
		}
		return false;
	};

	bool isInserted(int key, int value) {
		if (!isInserted(key)) return false;
		auto list = entries[key];
		return binary_search(list.begin(), list.end(), value);
	};

	vector<int> operator[] (int index) {
		//BT_PROFILE("MultiMap::operator[]");
		return entries[index];
	}
};


class ParticleSystem {
	
private:

	Vector3 ParticleSystemPosition;

	int maxNumParticle;
	int numParticle;

	float rho0;
	int resolutionParticles, numOfParticlesInKernel;
	Vector3 screenDomain, domain, scalingFactor;

	float constantOfPoly6, sigma, k, mu;
	float constantOfGradientSpiky;
	float constantOfLaplacianViscosity;
	float induceFactor;

	float maxSpeed;

	Vector3 gravity;

	int numXgrid, numYgrid, numZgrid;

	float timeStepSize;
	float maxSubTimeStepSize;

	vector<Vector3> pos;
	vector<Vector3> vel;
	vector<Vector3> mm;
	vector<float> rho;
	vector<float> pressure;
	Vector3 externalForce;

	vector<Vector3> worldPos;

	vector<Vector3> prs;
	vector<Vector3> vis;
	vector<Vector3> mfs;
	vector<Vector3> frs;
	vector<Vector3> tof;

	MultiMap neighborMap;
	MultiMap gridMap;
	
	vector<int> contactID;

	float m, h, h2, dt;

	vector<Vector3> aabbs[2];
	int numObjects;

	btDynamicsWorld* m_dynamicsWorld;

	Vector3 fixAxis;

	float maxStaticFriction2, dynamicStatic;

	float radius;	// minimum distacne between particle
	float radius_;
	
public:

	ParticleSystem(int _numParticle = 200) : domain(Vector3(1.0, 1.0, 1.0)), screenDomain(Vector3(8, 8, 8)),
		rho0(1000.0f), k(1000.0f), mu(10.0f), sigma(0.997f), maxSpeed(5.0), maxStaticFriction2(1500.0), dynamicStatic(200.0),
		gravity(Vector3(0, -9.8, 0)), numOfParticlesInKernel(10), maxNumParticle(1000),
		resolutionParticles(10000), maxSubTimeStepSize(0.0001), timeStepSize(0.001), fixAxis(1, 0, 1)
	{
		numParticle = _numParticle;
		scalingFactor = screenDomain / domain;
		m = rho0/(float)resolutionParticles;
		h2 = numOfParticlesInKernel*m / ((4.0f/3.0f)*rho0*PI);
		h = sqrt(h2);
		
		radius = h / 8.0;
		radius_ = radius * 0.95;

		constantOfPoly6 = 4.0f/(float)(PI*pow(h, 8));
		constantOfGradientSpiky = -30.0f/(float)(PI*pow(h, 5));
		constantOfLaplacianViscosity = 40.0f/(float)(PI*pow(h, 5));

		numXgrid = (int) (domain.x() / h) +1;
		numYgrid = (int) (domain.y() / h) +1;
		numZgrid = (int) (domain.z() / h) +1;
		
		induceFactor = 3 * Magnet::inverseU0;

		ParticleSeeding();

		externalForce = gravity;

	};
	~ParticleSystem() {
		neighborMap.Clear();
		gridMap.Clear();
	};

	void Simulation();
	void Simulation(const float timeStep) { timeStepSize = timeStep; Simulation(); };
	void SetDynamicsWorld(btDynamicsWorld* dw) { m_dynamicsWorld = dw; };

	void renderParticle();
	void renderParticle(bool);

private:

	void ParticleSeeding();

	void Advection();
	void Neighboring();
	void ComputeRho();
	void ComputeForceAndVelocity();

	void ComputeMoment();

	void SpatialGrid(Vector3& pos, int index);

	void GetAabbFromWorld();
	void DetectCollision(Vector3& prePos, Vector3& pos, Vector3& vel);
	void DetectCollisionWithParticle();

	int GridIndex(const Vector3& v);
	int GridIndex(int x, int z);
	int GridIndex(int x, int y, int z);

	const float B_spline_kernel(Vector3 xi, Vector3 xj, float h);
};


#endif